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ABSTRACT  Motivated  by  photodissociation  experiments  in  which  non-RRKM  nano¬ 
second  lifetimes  of  the  ethyl  radical  were  reported,  we  have  performed  a  classical  trajectory 
study  of  the  dissociation  and  isomerization  of  C2H5  over  the  energy  range  100—150  kcal/ 
mol.  We  used  a  customized  version  of  the  AIREBO  semiempirical  potential  (Stuart,  S.  J.;  et 
al.  /.  Chem.  Phys.  2000,  112,  6472— 6486)  to  more  accurately  describe  the  gas-phase 
decomposition  of  C2H5.  This  study  constitutes  one  of  the  first  gas-phase  applications  of  this 
potential  form.  At  each  energy,  10  000  trajectories  were  run  and  all  underwent  dissociation  in 
less  than  100  ps.  The  calculated  dissociation  rate  constants  are  consistent  with  RRKM 
models;  no  evidence  was  found  for  nanosecond  lifetimes.  An  analytic  kinetics  model  of 
isomerization/dissociation  competition  was  developed  that  incorporated  incomplete  mode 
mixing  through  a  postulated  divided  phase  space.  The  fits  of  the  model  to  the  trajectory  data 
are  good  and  represent  the  trajectory  results  in  detail  through  repeated  isomerizations  at  all 
energies.  The  model  correctly  displays  single  exponential  decay  at  lower  energies,  but  at 
higher  energies,  multiexponential  decay  due  to  incomplete  mode  mixing  becomes  more 

apparent.  At  both  ends  of  the  energy  range,  we  carried  out  similar  trajectory  studies  on  CD2CH3  to  examine  isotopic  scrambling. 
The  results  largely  support  the  assumption  that  a  H  or  a  D  atom  is  equally  likely  to  dissociate  from  the  mixed-isotope  methyl  end 
of  the  molecule.  The  calculated  fraction  of  products  that  have  the  D  atom  dissociation  is  ~20%,  twice  the  experimental  value 
available  at  one  energy  within  our  range.  The  calculated  degree  of  isotopic  scrambling  is  non-monotonic  with  respect  to  energy 
due  to  a  non-monotonic  ratio  of  the  isomerization  to  dissociation  rate  constants. 


1.  INTRODUCTION 

The  dissociation  of  the  ethyl  radical  (C2H5)  plays  an  important 
role  in  combustion  chemistry.  Because  the  reverse  reactions 
constitute  the  addition  of  a  hydrogen  atom  to  a  stable  molecule 
(as  opposed  to  radical— radical  recombination),  the  bond 
energies  are  relatively  small.  Consequently,  this  radical  tends 
to  dissociate  rapidly  in  the  primary  reaction  zones  of  premixed 
flames.  The  hydrogen  atom  thus  produced  acts  to  promote 
chain  branching  through  the  H  +  02  ^  OH  +  O  reaction. 
Thus,  the  rate  coefficients  for  this  dissociation,  along  with  those 
of  several  other  small  radicals,  are  important  factors  in 
determining  such  global  combustion  properties  as  laminar 
flame  speeds,  induction  times,  and  flammability,  extinction,  and 
ignition  limits.1,2 

Because  of  its  importance  in  combustion,  the  thermal 
kinetics  of  ethyl  radical  dissociation  has  been  studied  for 
many  years,3  and  the  most  recent  reviews4-6  show  a  consensus 
for  the  value  of  the  thermal  rate  constant  as  functions  of 
temperature  and  pressure.  This  consensus  includes  strong 
agreement  between  experiment  and  theory7  for  the  dissociation 
kinetics. 

In  addition  to  thermal  kinetics  experiments,  there  have  been 
several  photoinduced  dissociation  experiments  that  access 
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either  the  lower  3s  Rydberg  band8  11  in  C2H5  or  the  higher 
3p  band,12  or  both.13  These  experiments  create  ethyl  radicals  by 
photolysis  or  pyrolysis  from  some  suitable  precursor,  activate 
the  radical  by  photoabsorption  of  260  to  193  nm  light,  and  then 
photodetect  the  H  atoms  in  various  ways.  Before  describing  the 
experimental  results,  a  brief  review  of  the  experiments  is  useful. 
There  are  four  potential  sources  of  C2H5  in  these  experiments. 
The  photoabsorption8'9'12  of  the  ethyl  radical  places  it  in  a 
Rydberg  state  from  which  a  H  atom  can  directly  dissociate 
adiabatically,  leaving  behind  an  excited  C2H4.  The  H  atom  can 
also  migrate  between  the  carbon  atoms,  passing  through  a 
conical  intersection  to  the  ground  electronic  state.  (Theoretical 
calculations  of  reaction  paths  on  the  electronically  excited-state 
potential  energy  surfaces  (PESs)  of  C2H514  predict  that  the 
nonclassical  bridge  structure  is  a  minimum  on  the  3s  PES  and  a 
conical  intersection  to  the  ground  state.  The  bridge  structure  is 
also  a  minimum  in  the  3p  state  but  does  not  intersect  with  the 
ground  state.)  Once  in  the  ground  state,  C2HS  can  dissociate  to 
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a  H  atom  and  ground-state  C2H4,  either  directly  or  indirectly 
after  temporary  residence  in  the  well  on  the  ground-state  PES. 
A  fourth  source8  of  H  atoms  may  come  from  the  precursor  that 
is  a  monosubstituted  ethane  where  the  substituent  R  is  weakly 
bonded.  The  radical  R  can  then  attack  the  ethyl  radical,  form 
RH,  and  subsequently  be  photodissociated  along  with  the  ethyl 
radical.  Alternatively,  R  itself  can  dissociate  and  its  products  can 
be  subsequently  photodissociated.  While  the  characteristics  of 
the  dissociated  H  atoms  from  all  four  sources  are  of  intrinsic 
dynamical  interest,  the  only  source  of  particular  relevance  to 
combustion  is  the  unimolecular  dissociation  on  the  ground 
electronic  state  surface. 

Only  the  experiments  on  the  3s  Rydberg  band  of  the  ethyl 
radical  have  produced  results  characteristic  of  a  slow 
unimolecular  dissociation  of  the  radical.  The  experimental 
results  (including  a  lack  of  consensus  on  certain  points)  may  be 
summarized  as  follows: 

(1)  ~80%  of  photoexcited  C2H5  dissociates  by  the 
unimolecular  process/  However,  Zierhut  et  al.10  found  only  a 
fast  (20  fs)  decay  of  C2H5,  which  indicates  direct  dissociation; 
they  did  not  observe  a  slower  dissociation  that  could  involve  an 
intersystem  crossing.  This  interpretation  assumes  that  vibra- 
tionally  excited  C2H5  on  the  ground  state  surface  lives  long 
enough  to  be  photoionized  and  subsequently  detected  by  mass 
spectroscopy.  Recently,  Steinbauer  et  al.11  confirmed  the 
existence  of  both  a  unimolecular  and  direct  dissociation 
channel. 

(2)  The  differential  cross  section  for  that  unimolecular  decay 
is  isotropic,9  consistent  with  an  excited  complex  that  lasts  at 
least  a  rotational  period  (~1  ps). 

(3)  The  fraction  of  energy  released  into  translation  (fr)  is 
unexpectedly  high  from  earlier  work:  0.27  from  ref  8  and  0.35 
from  ref  9.  The  most  recent  study1 1  using  velocity  map  imaging 
has  (fr)  =  0.19  for  H  atoms  produced  by  the  slow  unimolecular 
process.  This  study  suggests  several  experimental  issues  that 
could  account  for  the  higher  values  of  earlier  studies.  The 
earliest  study8  characterized  the  translational  energy  distribu¬ 
tion  as  Maxwellian.  Subsequent  studies  have  shown  a  bimodal 
distribution  of  translational  energy  of  the  H  atoms  due  to  direct 
and  indirect  dissociation.  Decomposition  of  the  distribution 
shows  a  low  energy  distribution  that  is  Maxwell-like  but  is  well 
represented  by  a  more  involved  distribution.11 

(4)  The  unimolecular  dissociation  rate  is  multiexponential 
with  a  slow  component  of  106— 107  s_1  depending  on  the 
photon  energy.8'11  Rice  — Ramsperger  — Kassel— Marcus 
(RRKM)  calculations8  for  the  energy  range  of  the  experiment 
predict  a  rate  of  1012  s-1.  The  5  orders  of  magnitude  difference 
in  the  statistical  theory  rate  and  the  measured  rate  is  quite 
surprising  and  of  course  suggests  nonstatistical  behavior  at  very 
high  energies.  This  slow  rate  is  derived  from  the  measured 
temporal  growth  and  decay  (due  to  diffusion  out  of  the  laser 
focus)  of  H  atoms  over  ~400  ns.  The  measurements  that 
produce  these  slow  rates  do  not  have  the  time  resolution  to 
observe  picosecond  decay  rates.  Thus,  these  measurements  can 
only  be  interpreted  to  mean  that  there  is  a  component  of 
unimolecular  decay  many  orders  of  magnitude  slower  than  that 
predicted  by  RRKM. 

Additional  experimental  results  for  isotopic  scrambling  are 
available.  C2HS  can  both  dissociate  and  isomerize  (i.e.,  undergo 
a  1—2  hydrogen  shift): 

H2CCH3  -+  H2CCH2  +  H  (la) 


H2CCH3  ^  H3CCH2  (lb) 

The  isomerization  process  (eq  lb)  is  only  detectable  through 
isotopic  scrambling: 

D2CCH3  ->•  HD,CCH2  HDCCH2  +  D  (2) 

The  earliest  photodissociation  experiment13  that  accesses  the 
3s  Rydberg  band  found  some  isotopic  scrambling  by  detecting 
for  CD2CH3  one  D  atom  for  every  three  H  atoms.  However, 
this  early  experiment  used  193  nm  light  to  photodissociate  the 
precursor  and  that  could  contaminate  the  measured  isotopic 
ratio.  Two  subsequent  photodissociation8’9  experiments  that 
access  this  band  did  not  detect  any  isotopic  scrambling; 
however,  the  most  recent  experiment11  did  detect  isotopic 
scrambling  at  a  level  of  no  more  than  10%.  All  experiments12,13 
that  directly  access  the  3p  Rydberg  band  show  isotopic 
scrambling  that  is  attributed  to  isomerization  on  the  excited 
state  surface.12  In  a  thermal  dissociation  experiment15  of 
CD2CH3  and  CH2CD3  over  the  temperature  range  from  ~700 
to  ~1200  K,  only  the  majority  isotope  was  detected  as  the 
atomic  dissociation  product.  The  earliest  but  indirect 
scrambling  study16  involved  D3CCH3  pyrolysis  in  which  it 
was  argued  that  product  distributions  require  that  isomerization 
be  ~20%  of  the  dissociation  rate.  It  would  appear  that  the 
current  consensus  is  that,  on  the  ground  state  PES,  isomer¬ 
ization  is  at  best  a  minor  process  relative  to  dissociation  at  both 
thermal  and  photoexcitation  levels  of  energy 

The  history  of  theoretical  studies  '  relevant  to  the 
dynamics  of  C2H5  on  the  ground  state  PES  mostly  precedes 
all  the  photodissociation  methods  by  several  decades.  The 
earliest  trajectory  studies  by  Hase’s  group  using  analytical17'19,20 
PESs  based  on  ab  initio  electronic  structure  calculations  found 
at  higher  energies  (>100  kcal/mol)  that  unimolecular  decay  is 
described  by  a  single  rate  constant  well  represented  by 
statistical  theories  and  of  the  order  of  1012  s_1.  Their  computed 
product  translational  energy  distribution  has  the  average  value 
(fr)  =  0.18.  Isomerization  was  rarely  observed  in  the 
trajectories,  not  surprising  since  the  isomerization  barrier  is 
~20  kcal/mol  higher  than  the  dissociation  barrier  on  their 
PES17  and  isomerization  was  restricted  to  only  one  of  the  H 
atoms.17  Later,  electronic  structure  theory  calculations24 
predicted  this  barrier  to  be  only  5.0  kcal/mol  above  the 
dissociation  barrier.  Hase  and  co-workers21  subsequently  used 
high  quality  ab  initio  results  and  transition  state  theory 
calculations  to  obtain  good  agreement  with  measured  thermal 
rate  constants  for  dissociation  and  addition.  This  series  of 
calculations  are  either  consistent  or  unrelated  to  the  photo¬ 
dissociation  experimental  results  with  the  exception  of 
providing  no  evidence  for  a  very  slow  component  to  the 
unimolecular  decay. 

After  the  measured  results  of  Gilbert  et  al.8  of  the  slow  decay, 
three  different  direct  dynamics  trajectory  studies25  27  have 
been  carried  out  by  the  Chen  group  to  understand  the 
implications  of  the  measurements.  Two  of  those  studies25'26 
were  carried  out  at  the  HCTH147(3)6-31+G**/6-31G**  level 
of  theory  (see  the  Glossary  for  details  of  this  functional).  Those 
two  studies  found  biexponential  decay  that  they  associated  with 
trapped,  quasi-periodic  trajectories  (slow  rate)  in  a  sea  of 
chaotic  trajectories  (fast  rate).  A  wavelet-based  analysis28  was 
used  to  identify  specific  ethyl  vibrational  modes  involved  in  the 
quasi-periodic  motion.  The  third  study27  involved  surface¬ 
hopping  trajectories  propagated  on  the  multiple  electronic 
states  involved  in  the  photoabsorption  process.  State-averaged 
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multiconfiguration  self-consistent  field  (SA-MCSCF)  and 
multireference  configuration  interaction  with  singles  and 
doubles  (MR-CISD)  methods  had  to  be  used  to  provide  a 
balanced  description  of  both  ground  and  excited  electronic 
state  PESs.  This  study  found  that  only  13—25%  of  the  originally 
excited  ethyl  radicals  underwent  indirect  unimolecular  dis¬ 
sociation  on  the  ground  or  excited  PESs.  The  study  did  find 
that  those  trajectories  that  underwent  slow  unimolecular  decay 
on  the  ground  state  surface  had  a  largely  isotropic  distribution, 
had  a  (fT)  of  0.24,  had  a  single  unimolecular  dissociation  rate  of 
~0.5  X  1012  s“\  and  isomerized  on  average  ~1.25  times. 
Although  each  trajectory  typically  isomerized  at  least  once,  the 
dissociated  atom  always  started  on  the  methyl  end  of  the 
radical;  in  other  words,  no  scrambling  was  observed  even 
though  isomerization  occurred.  This  study  is  largely  consistent 
with  all  the  experimental  measurements  of  unimolecular  decay 
except  for  no  discovery  of  a  slow  decay  channel.  The  offered 
explanation  is  that  slow  decay  is  produced  by  trajectories 
trapped  by  the  near  resonance  of  vibrational  frequencies.  Such 
resonances  can  be  strongly  affected  by  curvature  and 
anharmonicity  changes  produced  by  different  electronic 
structure  methods.  Whether  any  of  the  three  direct  dynamics 
trajectory  studies  carried  out  by  the  Chen  group  are  definitive 
enough  to  reliably  describe  resonances  is  open  to  question. 

Of  more  practical  concern,  all  three  of  the  Chen  group 
studies  involved  only  hundreds  of  trajectories  propagated  on 
the  ground  state  PES  for  at  most  7  ps.  (The  surface  hopping 
study  comprised  5000  trajectories,  but  unimolecular  dynamics 
on  the  ground  state  surface  was  only  accessed  by  ~15%  of  the 
trajectories.  Relevant  to  the  previous  discussion  of  the  four 
sources  of  H  atoms  in  the  experiment,  this  simulation  suggests 
that  the  unimolecular  dissociation  is  a  minor  channel.)  Since 
the  experimental  observation  is  a  100  ns  lifetime  for  the  slow 
dissociation  channel,  the  simulations  are  a  factor  of  10  000  too 
short  to  confirm  experiment.  Since  only  a  small  fraction  of 
trajectories  displayed  the  slow  process,  uncertainty  from 
statistical  sampling  alone  barely  allows  the  slow  process  to  be 
distinguished  from  the  fast  process  over  such  a  short  time 
interval. 

As  discussed  above,  because  of  the  inherent  expense  of  direct 
dynamics,  all  of  the  theoretical  studies  designed  to  better 
understand  the  photoinduced  ethyl  radical  dissociation  experi¬ 
ments  were  based  on  too  few  trajectories  of  too  short  duration. 
There  are  only  two  ways  to  lift  these  restrictions.  First,  one  can 
fit  electronic  structure  calculations  on  the  PES  to  a  computa¬ 
tionally  inexpensive  but  accurate  functional  form  and  run  more 
and  longer  trajectories  with  that  form.  Second,  one  can  use  an 
inexpensive  semiempirical  PES  for  the  trajectory  runs.  We  are 
interested  in  both  approaches,  and  this  article  reports  our 
results  for  the  second  approach. 

Even  with  inexpensive  PESs,  running  many  trajectories  for 
hundreds  of  nanoseconds  is  a  formidable  endeavor.  In  the 
present  case,  the  dynamics  is  really  the  competition  between 
two  elementary  processes  (dissociation  and  isomerization)  that 
can  be  readily  defined  and  identified  in  the  trajectories.  This 
suggests  that  long-time  dynamics  could  be  followed  by  a 
kinetics  model  in  which  trajectories  determine  the  rates  of  the 
elementary  processes  and  the  solution  to  the  differential 
equations  of  the  model  determines  the  long-time  behavior. 
There  is  another  advantage  to  using  a  kinetics  model.  If  the 
dynamics  is  fully  statistical,  then  a  kinetics  model  will  display 
only  one  rate  for  product  formation.  However,  if  the  dynamics 
is  not  fully  statistical,  then  the  rates  of  other  processes,  such  as 


intramolecular  vibrational  redistribution  (IVR)  of  energy,  can 
influence  the  effective  dissociation  and/or  isomerization  rates. 
In  such  circumstances,  kinetics  models  will  display  multi¬ 
exponential  behavior.  In  other  words,  the  competition  between 
dissociation  and  isomerization,  as  well  as  dissociation  and/or 
isomerization  by  itself,  can  display  how  divided  phase  space  is. 

The  rest  of  this  paper  describes  a  trajectory  study  in  which 
10  000  trajectories  at  several  energies  were  followed  up  to  100 
ps  on  a  semiempirical  PES.  The  trajectories  where  converted 
into  “populations”  of  unisomerized  C2H5,  once  isomerized 
C2H5,  twice  isomerized  C2H5,  etc.  Various  kinetics  models  are 
fit  to  the  temporal  behavior  of  the  normalized  populations. 
Both  the  trajectories  themselves  and  the  kinetics  models  are 
used  to  interpret  the  unimolecular  dynamics.  The  rest  of  this 
article  is  organized  as  follows.  In  section  2,  we  describe  the 
semiempirical  PES  and  the  computational  details  of  the 
trajectories.  In  section  3,  we  present  and  discuss  our  trajectory 
results  for  six  total  energies  that  bracket  the  photon  energies 
used  in  the  experimental  photodissociation  studies  of  C2H5.  In 
section  4,  we  derive  various  kinetics  models  that  are  used  to 
analyze  the  trajectories.  Section  5  includes  a  summary  and  a 
brief  discussion  of  future  work. 

2.  COMPUTATIONAL  METHODS 

2.1.  PES  and  Level  of  Theory.  There  have  been  five 
previous21,25’”7  PESs  derived  from  ab  initio  or  semiempirical 
electronic  structure  calculations.  The  five  sets  of  calculated 
energetics  on  the  ground  state  PES  are  listed  in  Table  1,  and 


Table  1.  C2HS  Energies  in  kcal/mol  Calculated  by  the  Listed 
Electronic  Structure  and  Semi-Empirical  Methods  for  the 
Isomerization  and  Dissociation  Saddle  Points,  and  the 
Dissociation  Products  (the  Zero  of  Energy  is  the  Energy  of 
the  Equilibrium  Position  of  C2H5;  the  Bottom  Row  Gives 
the  Experimental  Dissociation  Energy) 


H 

A 

H2C-CH2 

C2H4—H 

C2H4('S)+H 

Present  Work 

CBS-QB3 

42.7 

40. 5a 

39.7 

M05-2X/CBSB7 

43.2 

45.7 

41.6 

AIREBO 

50.9 

58.1 

49.1 

AIREBO-SRP 

42.8 

40.3 

39.7 

Past  Work 

HCTH147@6-3  l+G**/6-3  lG**b 

45.9 

42.1 

CCSD(T)/CBSb 

42.4 

40.1 

QCISD(T)/6-3 1 1+G(2df,p)c 

41.6 

38.1 

MR-C 1 SD/6-3 1  ++G(d,p)° 

40.7 

35.9 

SA-MCSCF/6-3 1++G(d,p)d 

32.9 

25.8 

Experiment6 

39. 7f 

“See  ref  39.  ^References  25  and  31.  Reference  21.  ^Reference  27. 
^Reference  29.  ^Reference  30. 


four  of  the  sets  of  energetics  are  displayed  in  Figure  1.  The 
energy  scale  of  the  figure  is  expanded  to  include  the  energy 
range  of  the  photodissociation  experiments  of  Gilbert  et  al.,8 
Amaral  et  al./  and  Steinbauer  et  al.1  (all  in  the  darker  band) 
and  those  of  others1"'13  (the  lighter  band).  Both  the  table  and 
the  figure  show  only  classical  results  without  zero  point 
corrections.  (Since  quasi-classical  trajectories  are  used  to  follow 
the  dynamics,  initial  zero  point  energy  is  not  retained  by  the 
trajectories  during  the  unimolecular  dissociation  event. 
Consequently,  displaying  energetics  without  zero  point 
corrections  is  most  relevant  to  the  trajectory  dynamics 
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Figure  1.  Calculated  energetics  diagram  of  C2H5  including  both  the 
dissociation  and  isomerization  transition  state  and  the  dissociation 
asymptote.  The  different  methods  of  calculation  are  described  in  the 
text.  The  one  experimental  result  in  the  figure  is  labeled  ATcT  and  is 
described  in  the  text. 


described  in  section  3.)  The  most  recent  experimental 
determination29  of  D0  for  C2H5  — »  C2H4  +  H  is  34.9  kcal/ 
mol.  For  Figure  1,  De  was  obtained  from  D0  via  high  level 
anharmonic  vibrational  frequency  calculations30  (as  the 
experimentally  measured  vibrational  frequencies  for  C2H5  are 
incomplete)  to  obtain  the  39.7  ±  0.1  kcal/mol  value  listed  in 
Table  1.  The  five  sets  of  energetics  in  the  literature  can  be 
placed  in  two  groups.  As  discussed  in  section  1,  one  group 
(MR-CISD  and  SA-MCSCF)  was  part  of  a  study27  of  coupled 
electronic  surfaces  where  the  emphasis  was  on  a  balanced 
description  of  two  PESs,  resulting  in  a  compromise  in  the 
accuracy  for  the  ground  state  PES  which  is  of  relevance  here. 

Their  dissociation  energy  is  lower  (and  in  the  case  of  SA- 
MCSCF  very  much  lower)  than  experiment.  As  mentioned  in 
section  1,  the  second  group  has  a  natural  progression:  the 
energetics  of  quadratic  configuration  interaction  with  single, 
double,  and  perturbative  triple  [QCISD(T)]  with  a  large  basis 
set,21  the  energetics  of  coupled  cluster  with  single,  double,  and 
perturbative  triple  [CCSD(T)]  with  a  complete  basis  set  (CBS) 
extrapolation,31  and  energetics  of  the  HCTH147  density 
functional  theory  (DFT)-based  semiempirical  molecular  orbital 
method25  with  specific  reaction  parameters  trained  to  the 
CCSD(T)/CBS  results.  The  most  accurate  calculation,  CCSD- 
(T)/CBS,  agrees  very  well  with  the  measured  dissociation 
energy.  As  discussed  in  section  1,  the  least  expensive  of  these 
five  methods  is  still  too  expensive  to  implement  in  an  ab  initio 
molecular  dynamics  calculation,  leading  to  only  a  few 
trajectories  followed  for  short  times.  None  of  these  calculations 
determined  the  isomerization  saddle  point  and  consequent 
energetics. 

In  order  to  determine  the  isomerization  barrier,  a  relatively 
high  level  complete  basis  set-quadratic  Becke3  (CBS-QB3)  was 
carried  out  for  all  of  the  species  in  Figure  1.  As  seen  in  Table  1 
and  Figure  1,  the  resulting  isomerization  barrier  is  2.5  kcal/mol 
higher  than  the  resulting  dissociation  barrier.  This  is  about  half 
the  difference  calculated  by  a  pioneering  but  less  accurate 
polarization  configuration  interaction  (POL-CI)  approach  30 
years  ago.24  At  the  energies  of  the  photodissociation  experi¬ 
ments,  isomerization  and  dissociation  will  be  competing 


processes.  The  CBS-QB3  compares  quite  well  with  the 
experimental  D0.  The  CBS-QB3  along  with  the  only  other 
CBS  calculation  [CCSD(T)]  have  smaller  barriers  for  addition 
of  H  to  C2H4  than  any  of  the  other  methods  in  Table  1  or 
Figure  1.  Small  barriers  have  been  found  to  be  consistent  with 
the  measured  addition  kinetics  rate  constants.7  The  CBS-QB3 
molecular  structure  and  frequencies  of  C2H5  at  equilibrium,  the 
dissociation  saddle  point,  and  the  isomerization  saddle  point 
can  be  found  in  Tables  2,  3,  and  4,  respectively.  CBS-QB3 


Table  2.  C2HS  Equilibrium  Geometry  and  Harmonic 
Frequencies  (cm-1)  Calculated  by  Selected  Methods0 


CBS-QB3 

M05-2X 

AIREBO 

AIREBO-SRP 

Geometry 

Rcc 

1.487 

1.485 

1.517 

1.517 

Rem 

1.103 

1.098 

1.090 

1.090 

Rem 

1.094 

1.090 

1.090 

1.090 

i'CHl 

1.083 

1.079 

1.089 

1.089 

112.0 

111.6 

110.9 

110.9 

4>i 

111.9 

111.7 

110.6 

110.6 

<pi 

120.9 

120.8 

121.8 

121.8 

108.1 

108.3 

107.6 

107.6 

e. 

117.7 

118.0 

116.2 

116.2 

«3 

172.8 

173.2 

174.9 

174.9 

Frequencies 

108 

112 

127 

114 

477 

486 

655 

589 

813 

820 

688 

619 

980 

995 

1086 

976 

1063 

1088 

1399 

1258 

1192 

1206 

1407 

1265 

1401 

1419 

1521 

1367 

1465 

1488 

1800 

1618 

1483 

1504 

1891 

1700 

1483 

1507 

1957 

1760 

2943 

3003 

2828 

2542 

3035 

3089 

2864 

2575 

3078 

3134 

2923 

2628 

3139 

3192 

2928 

2632 

3239 

3297 

2959 

2660 

“The  geometry  uses  the  coordinate  system  of  ref  21.  The  bond  lengths 
are  in  A,  and  the  angles  are  in  degrees. 


calculated  values  for  C2H5  at  equilibrium  are  quite  similar  to  the 
large  basis  set  QCISD(T)  results  of  ref  21  and  other  high  level 
calculations.-3  Experimental  values  for  geometries  are  not 
readily  available  and  measured  vibrational  frequencies  are 
incomplete  and  anharmonic,  making  comparisons  difficult. 

To  find  much  faster  approaches  for  calculating  the  PES  that 
would  enable  following  many  trajectories  for  long  periods  of 
time,  the  energetics  (both  dissociation  and  isomerization)  of 
one  semiempirical  molecular  mechanics  method  and  two 
semiempirical  molecular  orbital  (SEMO)  methods  different 
from  HCTH147  were  investigated.  The  two  new  SEMO 
methods  were  parametric  method  number  6  (PM6)32  and 
M05-2X/CBSB7  (see  the  Glossary  for  details  of  this  func¬ 
tional).33  Both  have  been  trained  to  experimental  data  drawn 
from  a  variety  of  chemical  species.  PM6  produced  the  largest 
dissociation  energy  (~12  kcal/mol  higher  than  experiment), 
the  lowest  isomerization  barrier,  and  no  saddle  point  for  H 
addition  to  C2H4.  It  will  not  be  discussed  further.  In  contrast,  as 
can  be  seen  in  Table  1  and  Figure  1,  M05-2X  with  a  relatively 
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Table  3.  C2H5  Dissociation  Saddle  Point  Geometry 
(Distance  in  A,  Angles  in  deg)  and  Harmonic  Frequencies 
(cm-1)  Calculated  by  Selected  Methods” 


CBS-QB3 

M05-2X 

AIREBO 

AIREBO-SRP 

Rcc 

1.335 

Geometry 

1.334 

1.376 

1.310 

RcH1 

2.238 

2.147 

1.625 

2.238 

RcH2 

1.084 

1.081 

1.106 

1.090 

RzH3 

1.085 

1.081 

1.090 

1.090 

<P  1 

104.8 

103.8 

107.0 

104.9 

(p2 

121.7 

121.4 

120.2 

124.5 

$  3 

121.7 

121.4 

123.5 

124.5 

o2 

116.5 

117.0 

115.4 

110.9 

0, 

116.7 

117.2 

112.9 

110.9 

«3 

179.0 

178.9 

179.2 

166.6 

-308b 

Frequencies 

-534 

-6702 

-367 

241 

291 

708 

27 

273 

326 

808 

272 

833 

846 

873 

653 

961 

987 

942 

843 

975 

1008 

1339 

902 

1052 

1064 

1399 

1399 

1238 

1258 

1530 

1410 

1364 

1380 

1530 

1543 

1470 

1491 

1673 

1556 

1655 

1682 

1743 

1618 

3129 

3180 

2799 

2566 

3141 

3191 

2847 

2594 

3202 

3260 

2873 

2642 

3230 

3285 

2948 

2643 

“The  geometry  uses  the  coordinate  system  of  ref  21.  The  bond  lengths 
are  in  A,  and  the  angles  are  in  degrees.  bSee  ref  39. 


Table  4.  C2H5  Isomerization  Saddle  Point  Geometry 
(Distance  in  A,  Angles  in  deg)  and  Harmonic  Frequencies 
(cm-1)  Calculated  by  Selected  Methods” 


CBS-QB3 

M05-2X 

AIREBO 

AIREBO-SRP 

Geometry 

Rcc 

1.492 

1.487 

1.650 

1.650 

^CHl 

1.297 

1.293 

1.300 

1.297 

^CH2 

1.081 

1.078 

1.079 

1.079 

^CH3 

1.081 

1.078 

1.079 

1.079 

54.9 

54.9 

50.6 

50.5 

02 

119.4 

119.2 

118.9 

119.0 

03 

119.4 

119.2 

118.9 

119.0 

02 

118.6 

119.1 

115.8 

115.8 

03 

118.7 

119.1 

115.8 

115.8 

«3 

196.0 

195.5 

155.6 

155.8 

Frequencies 

-1928 

-2082 

-1928 

549 

553 

319 

746 

763 

576 

749 

764 

726 

790 

799 

821 

1131 

1160 

897 

1190 

1206 

923 

1274 

1287 

1256 

1418 

1436 

1662 

1441 

1467 

1666 

2217 

2248 

2217 

3139 

3201 

2630 

3144 

3205 

2647 

3243 

3308 

2713 

3266 

3329 

2716 

“The  geometry  uses  the  coordinate  system  of  ref  21.  The  bond  lengths 
are  in  A,  and  the  angles  are  in  degrees. 


inexpensive  basis  set  produced  dissociation  energetics  essen¬ 
tially  identical  to  HCTH147  that  was  trained  specifically  to  the 
C2H5  system  (although  both  methods  are  a  few  kcal/mol  too 
high  relative  to  CBS-QB3).  M05-2X  produces  an  isomerization 
barrier  within  0.5  kcal/mol  of  that  for  CBS-QB3.  In  Tables 
2—4,  the  geometries  and  frequencies  of  M05-2X  are  contrasted 
with  those  of  CBS-QB3.  The  M05-2X  equilibrium  geometry 
and  frequencies  are  in  good  agreement  with  the  CBS-QB3 
results.  In  Tables  3  and  4,  M05-2X  remains  in  superior 
agreement  to  CBS-QB3  results.  Consequently,  of  the  SEMO 
methods,  M05-2X  is  the  method  of  choice.  Unfortunately, 
while  M05-2X  is  a  reasonable  representation  of  the  PES,  in 
practice,  direct  ab  initio  trajectory  calculations  with  M05-2X 
proved  too  slow  for  the  purposes  of  this  study. 

Since  semiempirical  molecular  mechanics  potentials  are 
evaluated  with  no  diagonalization  procedure  for  eigenvalues 
(as  is  required  in  SEMO  methods),  such  potentials  are 
computationally  inexpensive  to  evaluate  and  well  suited  to 
trajectory  studies.  For  this  reason,  we  have  tested  the  adaptive 
intermole cular  reactive  empirical  bond-order  (AIREBO)34’35 
semiempirical  molecular  mechanics  potential  using  the  most 
recent  code  supplied  by  the  developer  S.  J.  Stuart.36  For  our 
application,  our  AIREBO  code  produces  a  potential  energy  of 
the  form  described  in  the  original  AIREBO  reference  but  with 
an  improved  set  of  parameter  values  (and  with  extensions  that 
do  not  apply  to  the  C2H5  system).  AIREBO  is  designed  to  be 
able  to  describe  reactions  and  systems  far  from  equilibrium. 
However,  it  has  been  calibrated  for  condensed  media 


application.  The  C2H5  study  constitutes  the  first  application 
of  AIREBO  to  a  gas-phase  molecular  system.  The  AIREBO 
energetics  are  listed  in  Table  1  and  displayed  in  Figure  1,  while 
the  AIREBO  geometries  and  frequencies  are  listed  in  Tables  2, 
3,  and  4  for  the  equilibrium,  dissociation  saddle  point,  and 
isomerization  saddle  point. 

While  AIREBO  is  orders  of  magnitude  faster  to  evaluate  than 
ab  initio  or  SEMO  methods,  it  has  technical  problems  and 
accuracy  issues  that  make  it  inappropriate  for  direct  application 
to  reaction  1.  The  technical  problems  concern  two  different 
discontinuities  in  the  gradients  of  the  potential.  The  first 
discontinuity  arises  from  a  linear  interpolation  of  terms  in  the 
potential  from  a  precalculated  table  of  their  values  over  a  grid  of 
interatomic  distances.  As  the  interatomic  distance  at  the 
evaluation  point  passes  through  a  node  in  the  grid  of  distances, 
one  of  the  two  grid  nodes  that  bracket  the  evaluation  distance 
abruptly  changes,  resulting  in  a  discontinuity  in  slope.  Direct 
calculation  of  terms  in  the  potential  without  table  look-up 
eliminates  this  problem.  The  second  discontinuity  is  in  the 
behavior  of  AIREBO  switches  at  isomerization.  AIREBO  for 
our  application  constructs  the  potential  V  from  two  parts: 

V'W  =  VREB0(X)  +  VLj(X)  (3) 

where  x  is  the  geometry  at  which  the  potential  is  being 
evaluated,  Vrebo  *s  a  reactive  empirical  bond  order  term,  and 
VLj  is  a  Lennard-Jones  term  for  corrections  to  Vrebo  at  larger 
distances.  Incorporated  into  VLr  is  a  complicated  product  of 
three  switches.  For  symmetric  positions  of  the  core  H2CCH2 
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Figure  2.  (a)  The  AIREBO-2010  potential  energy  as  a  function  of  (y,z)  position  of  a  H  atom  in  the  mirror  symmetry  plane  with  respect  to  the 
H2CCH2  fragment.  The  y  axis  is  along  the  C— C  bond  with  the  origin  on  the  right  C  atom.  Blue  designates  energies  lower  than  the  H  +  C2H4 
asymptote  (in  brick  red).  Yellow  designates  the  first  ~10  kcal/mol  above  the  asymptote.  Green  designates  higher  energies,  while  gray  designates 
energies  higher  than  100  kcal/mol.  The  red  contour  at  58.1  kcal/mol  indicates  the  dissociation  saddle  point,  (b)  The  AIREBO-2010  potential  energy 
as  a  function  of  (y,z)  position  of  a  H  atom  in  the  mirror  symmetry  plane  with  respect  to  the  H2CCH2  fragment.  The  y  axis  is  along  the  C— C  bond 
with  the  origin  on  the  midpoint  of  the  C— C  bond.  Color  coding  as  in  part  a  except  the  two  red  contours  at  50.8  and  52.0  kcal/mol  indicate  the 
isomerization  saddle  point. 


during  isomerization,  this  triple  switch  does  not  enforce  a  zero 
gradient  for  the  isomerizing  H  atom’s  perpendicular  passage 
through  the  C— C  perpendicular  bisector  plane.  The  result  is  a 
gradient  discontinuity  in  that  C— C  plane.  How  to  remedy  this 
problem  requires  discussion  concerning  accuracy  issues. 

AIREBO  is  insufficiently  accurate  for  our  purpose,  as  can  be 
seen  from  Tables  1—4  and  Figure  2.  The  two  panels  in  Figure  2 
represent  a  two-dimensional  cut  defined  by  the  Cartesian 
position  of  one  H  atom  through  a  15  degree  of  freedom  PES. 
The  horizontal  axis  is  coincident  with  the  C— C  bond.  Both 
panels  are  for  geometries  constrained  so  that  the  plane  defined 
by  the  C— C  bond  and  the  one  H  atom  is  a  mirror  symmetry 
plane  for  the  remaining  H2CCH2  structure.  Six  degrees  of 
freedom  orient  two  H  atoms  relative  to  the  C— C  bond  in  the 
H2CCH2  structure.  If  those  two  atoms  are  on  the  same  side  of 
the  mirror  symmetry  plane,  then  the  H2CCH2  structure  is  fully 
defined  when  the  seventh  degree  of  freedom,  the  C— C  bond 
length,  is  specified.  Each  panel  represents  the  lowest  energy  of 
the  H  atom  with  the  H2CCH2  structure  optimally  relaxed.  The 
two  panels  differ  in  the  definition  of  the  origin  of  the  plot.  For 
Figure  2a,  the  origin  is  one  of  the  C  atoms  and  incorporates  the 
equilibrium  and  dissociation  saddle  point  of  C2H5.  For  Figure 
2b,  the  origin  is  the  midpoint  of  the  C— C  bond  and 
incorporates  the  isomerization  saddle  point.  While  both  panels 
display  a  continuous  potential  surface,  the  surface  does  have 
kinks  in  it  (i.e.,  discontinuities  in  gradient).  In  Figure  2b,  the 
kink  along  part  of  the  isomerization  dividing  surface  (i.e.,  y  =  0 
in  the  plot)  is  due  to  the  second  technical  problem  mentioned 
above.  (The  first  technical  problem  concerning  table  look-up 
was  corrected  before  making  these  plots.)  However,  there  are 
other  kinks,  especially  in  Figure  2a,  that  are  unconnected  to  the 
AIREBO  PES  but  instead  are  due  to  the  relaxation  of  the 
H2CCH2  structure  in  making  the  plot.  There  can  be  multiple 
minima  in  the  seven  degrees  of  freedom  that  define  this 
structure,  and  which  minimum  is  the  lowest  can  abruptly 
change  with  the  position  of  the  H  atom  in  the  plot. 


The  results  in  Tables  1—4  and  in  Figure  2  indicate  at  least 
four  significant  inaccurate  features  with  the  native  AIREBO 
PES: 

•  The  energy  of  the  C2H4  +  H  asymptote  is  ~10  kcal/mol 
higher  than  experiment. 

•  The  spatial  extent  of  the  bonded  region  of  the  PES  is  too 
condensed.  The  dissociation  saddle  point  occurs  at  Rchi 
distances  that  are  ~0.6  A  shorter  than  ab  initio 
calculations. 

•  The  dissociation  barrier  is  ~10  kcal/mol  higher  than  the 
already  too  high  C2H4  +  H  asymptote,  while  the 
isomerization  barrier  is  comparable  to  the  asymptote. 
This  will  result  in  a  bias  to  isomerization  inconsistent 
with  the  ab  initio  energetics  for  these  two  asymptotes. 

•  The  isomerization  region  contains  an  unrealistically  large 
repulsive  region  at  Rchi  distances  larger  than  the  saddle 
point  (the  gray  region  in  Figure  2b).  Within  this  region, 
the  peak  repulsion  exceeds  200  kcal/mol. 

While  these  features  are  important,  the  scale  of  inaccurate 
features  is  ~10  kcal/mol,  while  the  energy  scale  of  the 
photoexcitation  experiments  discussed  above  is  ~100  kcal/mol. 
Furthermore,  the  equilibrium  properties  compare  reasonably 
well  to  ab  initio  calculations.  The  equilibrium  geometry  has  as 
its  maximum  error  an  RqC  bond  length  that  is  0.03  A  too  long. 
The  harmonic  zero  point  energy  is  only  ~600  cm-1  too  high  or 
about  40  cm-1  per  mode.  The  frequencies  are  consistently  too 
high  below  2000  cm-1  and  consistently  too  low  above  2000 
cm-1. 

Taken  together,  the  results  above  suggest  that  the  general 
AIREBO  model  is  a  useful  starting  point  for  modifications 
specific  to  this  PES  that  could  reduce  the  major  inaccuracies 
and  address  the  discontinuity  in  forces  discussed  above.  In  the 
spirit  of  past  work  on  other  generic  functional  forms,37'38  the 
modified  PES  will  be  call  AIREBO-SRP  for  specific  reaction 
parameters.  AIREBO  has  tens  of  parameters,  multiple  switches 
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Figure  3.  (a)  The  AIREBO-SRP  potential  energy  as  a  function  of  (y,z)  position  of  a  H  atom  in  the  mirror  symmetry  plane  with  respect  to  the 
H2CCH2  fragment.  The  y  axis  is  along  the  C— C  bond  with  the  origin  on  the  right  C  atom.  Blue  designates  energies  lower  than  the  H  +  C2H4 
asymptote  (in  brick  red).  Yellow  designates  the  first  ~5  kcal/mol  above  the  asymptote.  Green  designates  higher  energies,  while  gray  designates 
energies  higher  than  100  kcal/mol.  Red  contours  at  40.0,  40.2,  40.4,  41.0,  42.0,  and  43.0  kcal/mol  indicate  the  dissociation  saddle  point,  (b)  The 
AIREBO-SRP  potential  energy  as  a  function  of  ( y,z )  position  of  a  H  atom  in  the  mirror  symmetry  plane  with  respect  to  the  H2CCH2  fragment.  The 
y  axis  is  along  the  C— C  bond  with  the  origin  on  the  midpoint  of  the  C— C  bond.  Color  coding  as  in  Figure  3a  except  the  red  contours  at  40,  41,  42, 
and  43  kcal/mol  indicate  the  isomerization  saddle  point. 


of  a  variety  of  forms,  and  pairwise  functions  whose  form  can  be 
readily  changed.  Deep  modifications  are  possible  but  would  be 
sensible  only  in  the  context  of  systematic  ab  initio  calculations 
to  provide  the  basis  for  such  changes.  As  the  first  use  of 
AIREBO  in  the  context  of  an  isolated  gas  phase  reaction,  our 
goal  will  be  more  modestly  directed  at  using  only  the  ab  initio 
information  in  Tables  1—4  and  consequently  only  limited 
modifications  will  be  explored. 

There  are  three  modifications  that  are  easy  to  implement. 
First,  the  table  lookup  process  is  replaced  by  direct  calculation, 
thus  eliminating  one  source  of  the  discontinuities  in  the  forces. 
Second,  the  calculated  potential  energy  as  measured  relative  to 
the  C2H5  equilibrium  position  can  be  scaled  so  as  to  recover  the 
correct  H  +  C2H4  experimental  asymptote  (see  Table  l).  The 
required  scale  factor  is  a  fraction  (0.808).  (Unfortunately,  it  acts 
to  reduce  the  equilibrium  frequencies  into  poorer  agreement 
with  ab  initio  calculations.)  Third,  the  Lennard-Jones  terms  in 
eq  3  can  be  eliminated.  This  has  the  double  effect  of  eliminating 
the  discontinuity  in  the  forces  with  passage  through  the 
isomerization  plane  and  drastically  reducing  the  unphysical 
repulsive  energies  in  the  isomerization  region,  as  seen  in  Figure 
2b.  In  native  AIREBO,  the  o  parameter  in  the  Lennard-Jones 
potential  for  the  C— H  interaction  is  set  to  ~2.8  A.  Although 
the  switch  for  the  Lennard-Jones  potential  is  complicated,  it  can 
begin  to  turn  on  at  C— FI  distances  more  than  1  A  less  than  a, 
deep  within  the  repulsive  r-12  region  of  the  potential.  For  this 
C2H5  application,  the  switch  is  not  subtle  enough  to  navigate 
through  that  repulsive  region  in  a  physically  realistic  way. 

These  three  modifications  result  in  a  PES  that  is  still  too 
spatially  constricted  with  dissociation  and  isomerization  barriers 
still  poorly  related  to  ab  initio  calculations.  To  proceed  further, 
the  Vrebq  in  eq  3  must  be  more  fully  described: 


^reboM  ~  2  2  S 
>■ 


min  > 

ry-rH 

max  min 

vr*  ~r>i  , 


i  + 


Q.„ 


(4) 


where  rj;-  refers  to  interatomic  distances  that  the  vector  x  is 
composed  of,  S  is  a  switch  from  1  to  0  over  a  finite  range 
determined  by  the  min  and  max  parameters  in  their  arguments, 
btj  are  complicated  functions  of  switches  and  parameters  that 
determine  the  geometry-dependent  bond  order  for  the  ij  bond 
and  include  dependences  on  dihedral  angles  (and  thus  internal 
rotations),  and  everything  else  are  AIREBO  parameters.  The 
spatial  constriction  of  the  native  AIREBO  PES  is  due  to  the  fact 
that  rCHmax  is  set  to  1.8  A.  Any  H  whose  C— H  distance  exceeds 
that  amount  ceases  to  interact  with  any  C  atom;  i.e.,  H  has 
dissociated  from  the  C2H4  fragment  and  is  now  in  the  H  + 
C2H4  asymptotic  region.  Simply  extending  rCHmax  releases  the 
spatial  constriction  at  the  price  of  much  more  inaccurate  saddle 
point  properties. 

Consistent  with  a  minimalist  strategy,  our  additional 
modifications  focus  on  only  S  for  C— H  interactions.  All  other 
parameters  and  switches  remain  the  same.  In  particular,  fr  can 
include  S  for  C— H  interactions  within  its  switches.  For  that 
application,  S  remains  unchanged.  Only  as  the  outer  switch  in 
eq  4  and  only  for  the  C— FI  interaction  will  S  be  modified. 
Furthermore,  rCHmm  =  1.3  A  will  not  be  changed,  meaning  the 
equilibrium  region  of  the  PES  will  not  be  altered  (other  than  by 
the  scaling  factor  discussed  above).  Only  the  outer  regions  of  S 
will  be  changed,  and  the  changes  will  occur  in  two  stages  as 
discussed  below.  For  both  stages,  the  CBS-QB3  saddle  point 
properties  in  Tables  3  and  4  will  guide  the  modifications. 

In  the  first  stage,  S  is  modified  to  represent  the  dissociation 
saddle  point  without  regard  for  the  implications  for  isomer- 
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ization.  From  Table  3,  the  saddle  point  location  has  -RCH  =  2.2 
A.  rCHmax  must  be  greater  than  2.2  A,  and  the  final  value 
selected  is  2.5  A  At  that  extended  length,  for  RCH  >  rCHmln, 
there  are  extensive  regions  of  the  C— C— H  angle  (p  ^  90° 
where  H  interacts  with  both  C’s,  including  within  the 
equilibrium  region.  This  greatly  complicates  modifying  S  to 
reproduce  CBS-QB3  saddle  point  properties  while  maintaining 
a  physically  reasonable  PES  for  (p  >  90°.  The  solution  to  this 
difficulty  is  to  make  S  a  function  of  RCH  and  (p .  If  i pn  is  the  angle 
for  H  bonded  to  the  near  carbon  and  (p{  is  the  angle  for  H 
bonded  to  the  far  carbon,  then  for  (pn  >  90°,  it  turns  out  that  (pf 
<  63°  for  IAh  <  2.5  A.  S(Rch,  <p  <  60°)  can  then  be  set  to  the 
native  AIREBO  S  with  rCHmax  =  1.8  A,  ensuring  that  the 
equilibrium  region  is  not  perturbed.  However,  for  S(IA H,  <p  > 
90°),  rCHmax  is  set  to  2.5  A  and  its  native  AIREBO  functional 
form  is  abandoned  in  favor  of  a  point  by  point  fit  over  an  (R^; 
(p)  grid  to  a  template  analytic  function  designed  to  (l) 
reproduce  the  location,  height,  and  two  frequencies  of  the  CBS- 
0133  dissociation  saddle  point  in  (Rqh>  <P~)  space,  (2)  be 
physically  reasonable  for  (p  >  90°,  and  (3)  join  smoothly  in 
value  and  slope  the  native  AIREBO  PES  at  Rcjj  =  rCHmln. 
Intermediate  angular  grid  lines  between  60°  <  (p  <  90°  are  filled 
out  by  a  cubic  switch  in  angle  that  match  in  value  and  angular 
slope  the  values  at  <p  =  60  and  90°  (note  the  angular  slope  at 
60°  is  by  construction  0.0).  The  final  two-dimensional  grid  of  S 
values  for  1.3  A  <  RCH  <  2.5  A  and  60°  <  (p  <  180°  is  fit  to  a 
two-dimensional  spline  with  the  condition  of  zero  slope  on  all 
four  boundaries  of  the  grid. 

There  are  three  technical  details  regarding  this  approach. 
First,  the  template  potential  is  only  2  degrees  of  freedom  with 
the  presumption  that  the  other  13  degrees  of  freedom  are  fully 
optimized.  Consequently,  iterative  rounds  of  geometry 
relaxation  are  required  to  get  the  final  values  of  S  on  the 
two-dimensional  grid.  Second,  to  converge  saddle  point 
properties,  the  grid  must  be  considerably  finer  in  the  region 
of  the  saddle  point  than  elsewhere.  Hence,  the  final  grid  is 
regular  but  not  uniform.  Third,  for  C2HS,  the  two  saddle  points 
are  in  the  mirror  symmetry  plane  that  lies  along  the  C— C  bond 
of  the  C2H4  fragment.  This  makes  a  (p  dependence  to  S 
particularly  useful.  However,  there  is  an  azimuthal  angle  that 
can  locate  H  out  of  the  mirror  plane.  The  absence  of  that 
dependence  means  that  the  out-of-plane  saddle  point  frequency 
of  either  saddle  point  will  only  be  accidentally  represented  in 
this  approach. 

Figure  3a  and  the  AIREBO-SRP  column  in  Tables  1—3 
display  the  result  of  modifications  to  native  AIREBO  up  to  this 
stage.  Table  1  shows  a  correct  experimental  asymptote  and  a 
correct39  CBS-QB3  energy  relative  to  that  asymptote.  As  a 
consequence  of  a  transcription  error  in  the  CBS-QB3 
dissociation  barrier  height,  the  matched  value  for  the 
AIREBO-SRP  PES  is  40.3  kcal/mol  rather  than  40.5  kcal/ 
mol,  as  now  listed  in  Table  1.  Table  2  shows  that  the 
equilibrium  location  is  unchanged  from  that  of  native  AIREBO, 
but  the  frequencies  are  all  reduced  because  of  the  scale  factor. 
Despite  this  reduction,  in  general,  the  frequencies  are  still  too 
high  at  the  low  end  and  too  low  at  the  high  end  relative  to 
CBS-QB3  values.  Table  3  shows  that  the  CBS-QB3  saddle 
point  values  for  R^i  and  <pi  values  are  exactly  reproduced  and 
the  other  geometric  parameters  are  not  far  off.  The  imaginary 
frequency  and  the  in-mirror-plane  frequency  of  the  dissociating 
H  are  also  in  agreement39  with  CBS-QB3.  Similarly  to  Table  1, 
as  a  consequence  of  a  transcription  error  in  the  CBS-QB3 
tunneling  frequency,  the  matched  value  for  the  AIREBO-SRP 


PES  is  —367  cm-1  rather  than  —308  cm-1,  as  now  listed  in 
Table  3.  The  out-of-mirror-plane  frequency,  while  smaller  than 
the  in-mirror-plane  frequency  in  agreement  with  CBS-QB3,  is 
an  order  of  magnitude  too  small  relative  to  the  CBS-QB3  value. 
This  is  a  consequence  of  the  absence  of  an  azimuthal  angular 
dependence  in  the  modified  switch.  The  remaining  frequencies 
are  too  high  at  the  low  end  and  too  low  at  the  high  end,  as  in 
Table  2  for  the  equilibrium  frequencies. 

The  last  stage  of  modifications  concerns  the  isomerization 
region.  The  angle  the  isomerizing  H  atom  makes  with  either 
end  of  the  C— C  axis  is  acute.  The  modification  strategy  above 
is  inapplicable  to  this  region  because  the  switch  must  be 
simultaneously  modified  for  two  C— H  interactions.  The  last 
stage  of  modifications  takes  advantage  of  the  fact  that  the 
isomerization  saddle  point  is  located  at  the  intersection  of  two 
symmetry  planes,  namely,  the  mirror  symmetry  plane  and  the 
C— C  perpendicular  bisector  plane.  Aong  the  line  of 
intersection,  the  two  C— H  interactions  are  identical.  Following 
the  process  used  in  the  dissociation  region,  the  switch  value 
needed  to  reproduce  a  template  potential  along  that  line  for  a 
fully  relaxed  C2H4  fragment  can  be  determined.  The  template 
potential  can  be  shaped  to  have  the  CBS-QB3  saddle  point 
properties  along  that  line,  to  merge  seamlessly  into  the  native 
AIREBO  PES  at  values  of  Rqj  that  are  interior  to  the  saddle 
point,  and  to  be  physically  reasonable  at  distances  much  larger 
than  the  saddle  point.  If  the  resulting  switch  value  is  call  Siso, 
then  S  must  transform  into  Siso  as  IRcih  —  Rc2hI  approaches 
zero  (where  Cl  and  C2  indicate  the  two  carbons).  This  is 
accomplished  with  a  new  switch  S'  operating  on  the  argument  q 
=  IRcih  —  Rc2hI/AR,  where  A R  is  the  range  over  which  S' 
applies.  If  Smodl  is  the  modified  switch  after  the  first  stage  of 
modifications  discussed  previously  and  if  Smod2  is  the  final 
switch  after  this  second  stage  of  modifications,  then 

Smod2  =  S  ( <2 ) SiSO  +  [l  —  S  (q)]Smodl  (5) 

for  the  Cl— H  and  C2— H  interactions.  Equation  5  requires  that 
S'(q  >  l)  =  0  and  S'(0)  =  1,  while  continuity  in  value  and 
derivative  require  that  the  derivative  S'(q)  =  0  for  q  =  0  or  1. 
This  can  be  accomplished  by  an  even  powered  fourth  order 
polynomial  in  q.  However,  in  this  scheme,  S'(q)  determines  the 
imaginary  frequency  at  the  saddle  point.  S'(q)  was  given  more 
flexibility  to  set  this  frequency  by  making  it  a  sixth  order 
polynomial  of  the  form: 

S'(q)  =  1  —  6q2  +  9q4  —  4q6  +  4wq2(l  —  q2)2  (6) 

where  w  is  a  parameter  that  alters  the  width  of  the  switch.  If 
0.75  <  w  <  1.5,  S'(q)  has  no  extraneous  extrema. 

Figure  3b  and  the  AIREBO-SRP  column  in  Tables  1  and  4 
display  the  result  of  modifications  to  native  AIREBO  that 
include  isomerization.  For  these  results,  AR  =  0.196  and  w  = 
0.95,  parameter  values  that  leave  unchanged  the  first  stage 
modifications  displayed  in  Figure  3a  and  Tables  1—3.  Table  1 
shows  a  correct  CBS-QB3  isomerization  saddle  point  energy 
relative  to  the  experimental  H  +  C2H4  asymptote.  Table  4 
shows  that  the  CBS-QB3  saddle  point  value  for  Rchi  is  exactly 
reproduced.  The  perpendicular  distance  of  the  saddle  point 
from  the  C— C  midpoint  (i.e.,  z  in  Figure  3b)  is  also  exactly 
reproduced  but  not  displayed  in  the  table.  Reproducing  those 
two  distances  results  in  a  (p1  value  that  is  somewhat  too  low, 
unlike  the  exact  reproduction  of  the  CBS-QB3  value  in  Table  3 
for  dissociation.  While  (p1  is  too  low,  Rcc  is  too  large  but  all  the 
other  distances  and  angles  are  as  approximately  correct  in  Table 
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4  for  isomerization  as  they  are  in  Table  3  for  dissociation.  The 
large  imaginary  frequency  and  the  large  C— H  stretch  frequency 
along  the  perpendicular  direction  (z  in  Figure  3b)  are  in  exact 
agreement  with  CBS-QB3.  Unlike  the  results  in  Table  3  for 
dissociation,  the  out-of-mirror-plane  frequency  for  the  isomer- 
izing  H  at  1256  cm-1  is  quite  close  to  the  CBS-QB3  value  of 
1274  cm-1.  Like  the  Table  3  results,  the  remaining  highest 
frequencies  are  too  low.  Unlike  the  Table  3  results,  the  lower 
frequencies  (<2000  cm-1)  are  a  mix  of  too  high  and  too  low. 

Overall,  the  AIREBO-SRP  developed  here  avoids  disconti¬ 
nuities  in  the  gradients,  displays  CBS-QB3/experimental 
energetics  for  dissociation/isomerization  processes,  has  the 
correct  general  spatial  extent,  and  approximates  the  CBS-QB3 
geometries  and  frequencies  of  the  equilibrium  and  two  saddle 
points.  For  these  three  structures,  the  root-mean-square  (rms) 
error  on  distances  relative  to  CBS-QB3  values  varies  from 
~0.015  A  for  equilibrium  and  the  dissociation  saddle  point  to 
0.08  A  for  the  isomerization  saddle  point  (due  to  the  0.16  A 
discrepancy  for  Rqc)-  The  rms  error  in  the  angles  varies  from 
~2.5  to  ~6.5°.  The  rms  error  in  frequencies  is  more  substantial, 
ranging  from  ~300  to  ~350  cm-1.  These  errors  tend  to  be  in 
the  same  direction  so  that  the  difference  from  CBS-QB3  for  the 
change  in  zero  point  energy  between  saddle  point  and  the 
equilibrium  geometry  is  ~90  cm-1  for  dissociation  and  ~550 
cm-1  for  isomerization  out  of  an  average  zero  point  energy  of 
~11000  cm-1.  The  AIREBO-SRP  PES  will  be  used  for  the 
trajectory  runs  whose  results  will  be  described  in  the  rest  of  this 
article.  As  a  measure  of  continuity  of  the  AIREBO-SRP 
gradients,  no  propagation  errors  are  encountered  for  traditional 
step  sizes  (~0.05  fs)  and  energy  conservation  errors  >0.7%  of 
the  total  energy  are  restricted  to  <1%  of  the  trajectories  run.  On 
the  native  AIREBO  surface,  energy  conservation  errors  of  >1% 
of  the  total  energy  occurred  for  ~60%  of  the  trajectories. 

2.2.  Initial  Conditions  and  End  of  Propagation 
Conditions.  The  trajectories  were  calculated  using  a  version 
of  GENDYN40  interfaced  with  AIREBO-SRP.  Microcanonical 
ensembles  of  10  000  trajectories  were  calculated  for  total 
energies  at  10  kcal/mol  intervals  over  the  range  100—150  kcal/ 
mol.  The  initial  conditions  were  prepared  by  starting  with  the 
ethyl  radical  at  the  AIREBO-SRP  equilibrium  geometry  (see 
the  last  column  in  Table  2).  A  Markov  walk41  of  2  million  steps 
was  carried  out  to  move  the  system  away  from  equilibrium  to  a 
random  phase  space  point.  Since  the  radical  is  symmetric,  the 
Markov  walk  was  restricted  on  only  one  of  the  double  wells  by 
placing  a  reflecting  barrier  at  the  saddle  point  for  the  migration 
of  a  H  atom  from  the  methyl  carbon  to  methylene  carbon.  At 
each  Markov  step,  all  of  the  Cartesian  coordinates  were 
changed  by  a  random  fraction  of  0.11  A  with  sampling  over  the 
range  of  2.0  A  for  each  coordinate;  these  parameters  gave 
acceptance/rejection  ratios  over  the  range  0.44—0.64.  The 
configuration  at  the  end  of  the  initial  randomization  walk  was 
used  to  generate  the  initial  conditions  of  the  first  trajectory,  and 
the  initial  conditions  for  subsequent  trajectories  were  obtained 
by  performing  10  000  Markov  steps  from  the  configuration 
used  to  generate  the  initial  conditions  of  the  previous  trajectory. 
The  initial  conditions  were  set  so  that  the  total  angular 
momentum  was  zero. 

The  trajectories  were  integrated  using  velocity  Verlet  with  a 
time  step  of  0.05  fs,  which  gave  an  integration  accuracy  that 
resulted  in  less  than  1%  of  the  trajectories  having  greater  than  1 
kcal/mol  drift  for  the  total  energy  of  150  kcal/mol.  The  forces 
were  obtained  by  central  differencing  with  the  distance 
increment  1.0  X  10-5  A.42  A  trajectory  was  ended  if  a  C— H 


bond  length  exceeded  2.5  A  (that  is,  dissociation  occurred)  or 
the  total  integration  time  exceeded  100  ps.  Isomerization  was 
defined  to  have  occurred  when  a  hydrogen  atom  crossed  the 
isomerization  plane  and  the  resulting  new  C— H  bond 
underwent  at  least  one  inner  and  one  outer  turning  point 
(i.e.,  at  least  one  vibration)  without  recrossing  the  isomer¬ 
ization  plane.  The  time  of  the  isomerization  was  taken  to  be  the 
time  at  which  the  barrier  was  crossed. 

3.  TRAJECTORY  RESULTS  AND  DISCUSSION 

As  described  above,  10  000  trajectories  were  followed  for  100 
ps  at  six  different  energies  that  span  the  energies  of  the 
photodissociation  experiments.  The  results  of  these  simulations 
are  organized  below  under  the  topics  of  dissociation/isomer- 
ization  dynamics,  product  translational  energy  distributions, 
and  isotopic  scrambling. 

3.1.  Dynamics  of  Dissociation  and  Isomerization.  In 

order  to  study  both  dissociation  and  isomerization,  the 
trajectories  can  be  represented  in  terms  of  the  time-dependent 
normalized  population  P,(t)  of  trajectories  that  at  time  t  have 
undergone  only  i  isomerizations  and  have  not  yet  dissociated. 
Initially,  no  trajectories  have  undergone  isomerization  so  P0(t) 
can  only  decay  with  time  by  isomerization  or  dissociation.  All 
other  Pj(t)  will  initially  grow  with  time  from  isomerization  of 
Pj_i(t )  and  decay  at  later  times  from  dissociation  or 
isomerization.  The  solid  lines  of  Figure  4  show  the  normalized 
populations  of  the  C2H5  isomers  that  have  had  0  to  3 
isomerizations,  for  the  energies  of  100,  130,  and  150  kcal/mol. 
(The  figure  also  shows  with  dash— dot  lines  kinetic  model  fits  to 
these  populations  that  will  be  discussed  in  the  next  section.) 
The  results  at  all  three  energies  are  relatively  similar  with 
primarily  only  the  time  scale  changing.  These  results  show  two 
important  features  of  the  dynamics.  First,  isomerization 
populations  are  on  the  order  of  20%  of  the  initial  population 
and  multiple  isomerizations  are  clearly  evident.  At  thermal 
energies  (at  least  60  kcal/mol  lower  than  any  of  these  results), 
isomerization  will  poorly  compete  with  dissociation  because,  as 
shown  in  Figure  1  and  Table  1,  the  isomerization  barrier  is  2.5 
kcal/mol  higher  than  dissociation. 

The  second  feature  of  these  results  is  the  absence  of  very 
long-lived  trajectories.  Even  at  the  lowest  energy  represented 
here,  all  the  trajectories  were  followed  until  they  dissociated.  All 
of  the  trajectories  dissociated  within  100  ps  integration  time. 
These  results  do  not  support  any  significant  amount  of  C2H5 
dissociating  with  a  rate  of  10  ns.  However,  this  does  not  mean 
that  the  isomer  populations  dissociate  with  a  single  exponential 
decay.  In  Figure  5,  the  P0(t *)  trajectory  results  are  plotted 
versus  an  energy-scaled  time  t*  along  with  optimal  single 
exponential  fits,  f*  is  f/(l  +  [0.1(150  —  E)]1'75)  and  is  solely 
designed  to  allow  P0(t )  for  E  =  100  and  130  kcal/mol  to  be 
plotted  with  sufficient  resolution  on  the  same  graph  as  P0(t )  for 
E  =  150  kcal/mol.  The  single  exponential  fit  is  for  a  unit 
prefactor  (to  preserve  a  normalized  population  in  the  fit)  and  is 
optimized  to  minimize  the  absolute  rms  error.  As  the  figure 
shows,  at  low  energies,  a  single  exponential  fit  is  a  good 
representation.  However,  at  the  higher  energies,  at  both  high 
and  low  values  of  P0,  there  are  noticeable  deviations  on  a 
maximum  scale  approaching  0.05. 

The  decay  of  P0(t)  is  due  to  both  isomerization  and 
dissociation.  Tests  at  E  =  150  kcal/mol  show  that  both  the 
dissociation  products  and  the  isomerization  products  from 
P0(t)  decay  have  a  comparable  deviation  from  a  single 
exponential  model.  Although  Figure  5  displays  behavior  for 
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Time  (ps) 

Figure  4.  Normalized  trajectory  populations  Pf(f)  (solid  lines)  and 
corresponding  kinetics  model  fits  for  i  =  0,  1,  2,  and  3  for  a  total  energy 
E  of  100  kcal/mol  (panel  a),  130  kcal/mol  (panel  b),  and  150  kcal/ 
mol  (panel  c).  See  text  for  details. 


Energy-Scaled  Time  (ps) 

Figure  5.  P0(t *)  versus  energy-scaled  time  f*  for  E  =  100,  130,  and 
ISO  kcal/mol.  See  text  for  details. 


only  P0,  the  overall  H  product  formation  rate  (which  is  what 
experiments  directly  measure)  also  shows  multiexponential 
behavior  at  high  energies  but  single  exponential  behavior  at  low 
energies.  However,  the  multiexponential  behavior  observed 
does  not  involve  multiple  rate  constants  that  differ  by  many 
orders  of  magnitude,  as  suggested  by  the  Chen8  and  Fischer11 
groups. 

All  the  trajectories  used  in  this  study  have  no  initial  rotational 
energy.  As  such,  no  meaningful  differential  scattering  cross 
section  can  be  calculated.  However,  as  seen  in  Figure  4,  at  all 
energies  studied,  the  C2HS  lifetime  exceeds  1  ps.  The  rotational 
period  of  C2HS  is  ~1  ps  or  less,  depending  on  the  degree  of 
rotational  excitation.  Thus,  the  expectation  is  that  the  angular 
scattering  would  be  nearly  isotropic,  in  agreement  with 
experimental  measurements. 

3.2.  Translational  Energy  Distributions.  Figure  6  shows 
the  product  translational  energy  distribution  and  the  corre- 


Trans.  Energy  (kcal/mol) 

Figure  6.  Product  translational  energy  distributions  from  trajectories 
and  from  the  corresponding  least-squares  fit  Maxwell  distribution. 

sponding  optimal  Maxwell  distribution  fit  for  the  highest  and 
lowest  energies  studied.  (fT),  the  fraction  of  product  energy 
released  into  translation,  is  listed  for  each  energy  in  Table  5 

Table  5.  (fT)  Calculated  by  Trajectories  and  Tmax  from  the 
Maxwell  Distribution  Fit  to  the  Product  Translational 
Energy  Distribution 


E  (kcal/mol) 

<fT> 

Fmax(K) 

100 

0.139 

3300 

110 

0.133 

3700 

120 

0.129 

4100 

130 

0.126 

4500 

140 

0.123 

4900 

150 

0.121 

5200 

along  with  the  temperature  derived  from  the  Maxwell 
distribution  fit.  As  the  energy  increases,  the  trajectory  and 
Maxwell  distributions  shift  to  higher  and  broader  peaks. 
However,  the  Maxwell— Boltzmann  distribution  is  not  peaked 
enough  to  well  represent  trajectory  results.  From  Table  5,  (/T) 
varies  by  ~0.02  from  0.14  at  100  kcal/mol  to  0.12  at  150  kcal/ 
mol.  Hase  et  al.19  calculated  a  value  of  0.18  at  100  kcal/mol  on 
a  different  PES.  The  most  recent  and  probably  most  reliable 
experimental  measurement11  is  0.19.  As  discussed  earlier,  the 
product  translational  energy  distribution  is  bimodal,  of  which 
only  the  lower  energy  component  is  relevant  to  unimolecular 
decay  on  the  ground  state  PES.  Consequently,  the  experimental 
value  can  be  complicated  by  the  decomposition.  However,  it  is 
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most  likely  the  case  that  the  AIREBO-SRP  PES  results  in  a 
somewhat  low  value  for  (fT). 

3.3.  Isotopic  Scrambling.  For  both  ends  of  our  energy 
range,  trajectory  studies  analogous  to  Figure  4a  and  c  were 
carried  out  for  CD2CH3.  The  fraction  Fu  of  products  that 
constitutes  H  +  CD2CH2  was  determined  and  is  reported  in 
Figure  7  as  the  red  symbols;  it  is  substantially  less  than  unity, 
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Figure  7.  The  fraction  of  dissociations  of  CD2CH3  that  produce  H 
atoms  as  a  function  of  energy.  Trajectory  calculations  for  CD2CH3 
(red  symbols),  trajectory  calculations  for  CH2CH3  with  interpolated 
results  from  CD2CH3  (blue  symbols),  and  kinetics  model  fits  to 
CH2CH3  with  statistical  isotopic  fractions  (black  line).  See  text  for 
details. 


about  0.8.  In  contrast  and  as  shown  in  the  figure,  the  most 
recently  measured1  1  fh  in  this  energy  range  is  greater  than  0.9. 
This  PES  appears  to  allow  too  much  isotopic  scrambling  from 
isomerization.  Figure  8  examines  Fh  for  each  isomer  population 


c 


Figure  8.  As  a  function  of  the  number  of  isomerizations,  the  fraction  of 
dissociations  of  CD2CH3  that  produce  H  atoms  from  a  statistical 
model  and  from  trajectories  at  E  =  100  and  150  kcal/mol.  Error  bars 
on  trajectory  results  represent  one  o.  See  text  for  details. 


for  the  low  (red  line)  and  high  (blue  line)  energy.  For  early 
isomer  populations,  the  high  and  low  energy  results  alternate  in 
synchronization  between  high  and  low  H  fractions.  For  the  later 
isomer  populations,  the  sampling  statistics,  represented  by  one 
<7  error  bars,  become  too  poor  to  distinguish  the  two 
distributions.  Intermediate  isomer  populations  do  seem  to 
favor  more  H  dissociation  at  lower  energies. 

The  alternation  in  Figure  8  can  be  easily  understood  from  a 
statistical  model  in  which  H  and  D  can  equally  well  isomerize 
or  dissociate  from  the  methyl  end  of  the  radical.  With  such  an 
assumption,  CD,CH,  can  only  dissociate  to  H  atoms  giving  an 


Fn  of  1.  However,  CD2CH3  can  isomerize  to  produce 
CH2CHD2  whose  FH  is  1/3.  CH2CHD2  can  isomerize  to 
produce  67%  of  the  time  CH2DCHD  (Fh  of  2/3)  and  33%  of 
the  time  CH3CD2  (Fh  of  l).  This  results  in  an  FH  of  7/9.  The 
generalization  of  this  process  produces  an  Fh  of  the  ith  isomer, 

Fn(i), 

Fh(i)  =  0.6  +  0.4(— 2/3)'  (7) 

which  is  the  black  line  in  Figure  8.  Dissociation  after  a  very 
large  number  of  isomerizations  gives  an  Fh0»1)  =  0.6  or  the 
ratio  of  H  atoms  to  all  hydrogenic  atoms  in  CD2CH3,  the 
expected  result  for  complete  isotopic  scrambling. 

The  trajectory  results  always  show  a  Fh(0  higher  or  at  best 
the  same  as  the  statistical  result  of  eq  7,  implying  that  H  and  D 
do  not  behave  equally  with  regard  to  isomerization  and 
dissociation.  H  atoms  are  consistently  more  prone  to 
dissociation  than  D  atoms,  and  that  propensity  might  increase 
as  the  energy  decreases.  While  that  is  the  conclusion  to  draw 
from  Figure  8,  this  conclusion  may  be  sensitive  to  the  definition 
of  isomerization.  As  discussed  earlier,  isomerization  here  occurs 
when  the  isomerizing  atom  crosses  the  isomerization  plane  and 
undergoes  an  outer  and  an  inner  turning  point.  By  this 
definition,  a  H  atom  that  crosses  the  isomerization  plane  and 
then  directly  (i.e.,  no  outer  turning  point)  dissociates  without 
ever  recrossing  the  isomerization  plane  would  be  said  to 
dissociate  without  isomerizing.  At  E  =  150  kcal/mol,  the  C2H5 
trajectories  of  Figure  4c  were  rerun  with  termination  upon 
either  dissociation  or  upon  an  isomerization  defined  solely  by 
crossing  the  isomerization  plane.  Comparison  with  the  original 
trajectory  run  indicates  that  ~20%  of  the  original  P0 
dissociations  crossed  the  isomerization  plane  at  least  once.  If 
a  substantial  number  of  those  trajectories  crossed  the  plane 
only  once  on  the  way  to  directly  dissociate  with  no  further 
turning  points,  then  an  isomerization  definition  based  solely  on 
crossing  the  isomerization  plane  would  result  in  the  hydrogenic 
atom  that  just  isomerized  having  a  greater  propensity  to 
dissociate  than  the  two  other  hydrogenic  atoms  already  resident 
on  that  side  of  the  radical.  This  would  violate  the  statistical 
assumption  underlying  eq  7  and  perhaps  quantitatively  change 
the  degree  of  alternation  seen  in  the  trajectory  results  in  Figure 
8. 

Applying  the  alternations  in  Figure  8  to  the  asymptotic 
product  population  for  each  isomer  from  the  isotopically  pure 
CH2CH3  trajectory  runs  are  shown  as  the  blue  triangles  in 
Figure  7.  In  this  application,  the  calculated  alternations  at  E  = 
100  and  150  kcal/mol  are  connected  by  a  linear  interpolation  in 
energy.  At  the  two  end  point  energies  where  no  interpolation 
error  is  involved,  the  conversion  of  CH2CH3  results  in  FH  for 
the  mixed  isotope  CD2CH3  having  less  than  a  0.02  error.  As 
seen  previously  in  Table  1,  the  best  available  electronic 
structure  calculations  have  an  isomerization  barrier  2.5  kcal/ 
mol  higher  than  the  dissociation  barrier.  This  implies  that,  as  E 
approaches  the  C— H  dissociation  energy  of  ~40  kcal/mol, 
trajectory  calculations  of  FH  will  approach  unity.  In  this  context, 
the  blue-triangle  trajectory  results  in  Figure  7  display  a  high- 
energy  minimum  or  at  least  a  plateau  in  the  energy  dependence 
of  Fh.  The  minimum  value  of  ~0.8  is  midway  between  the 
possible  extremes  of  1.0  (dissociation  dominates  isomerization) 
and  0.6  (isomerization  dominates  dissociation). 

As  discussed  earlier,  the  AIREBO-SRP  PES  was  only 
calibrated  to  a  variety  of  calculated  isomerization  and 
dissociation  saddle  point  properties.  The  total  energy  involved 
in  the  measured  Fh.  value  in  Figure  6  is  more  than  twice  the  C— 
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H  dissociation  energy  and  therefore  samples  regions  of  the  PES 
far  removed  from  saddle  points.  The  comparison  to  experiment 
in  Figure  7  indicates  that  the  region  of  the  isomerization  plane 
on  the  AIREBO-SRP  PES  is  too  “open”,  allowing  too  much 
isomerization  that  suppresses  substantially  below  the 
measured  value. 

4.  KINETICS  ANALYSIS  OF  TRAJECTORIES 

The  representation  of  the  trajectory  results  in  terms  of 
populations  (as  in  Figure  4)  calls  to  mind  coupled  kinetics 
models,  in  which  species  that  react  with  each  other  are  modeled 
by  first  order  kinetics  equations.  The  solutions  of  the  first  order 
kinetics  equations  give  the  species  concentrations  as  a  function 
of  time  that  parametrically  depend  on  the  reaction  rates.  The 
populations  displayed  in  Figure  4  have  the  property  of  a 
cascade  in  which  there  is  no  process  that  transforms  a 
population  P,(f)  into  a  “prior”  population  Pj<i(t).  The  kinetics 
of  cascades  can  often  be  represented  analytically,  leading  to 
representations  of  P,(t)  that  can  be  least-squares  fit  to  the 
trajectory  results  to  produce  optimal  values  of  the  rate 
parameters  in  the  kinetics  model.  This  section  describes  such 
a  kinetics  model.  In  addition  to  their  interpretive  value, 
validated  kinetics  models  can  be  extrapolated  to  higher  or  lower 
energies  and  to  longer  times  to  gain  a  more  comprehensive 
understanding  of  the  processes  studied  by  the  trajectories. 

A  successful  kinetics  model  must  be  able  to  describe 
dissociation  and  isomerization  processes  for  a  cascade  of 
isomers,  i.e.,  Figure  4.  For  quantitative  accuracy,  it  must  be  able 
to  describe  the  emergence  of  multiexponential  decay  with 
increasing  E,  i.e.,  Figure  5.  For  extrapolation  as  a  function  of 
energy,  each  optimized  parameter  value  at  the  six  different 
values  of  E  where  trajectories  were  calculated  must  have  a 
“smooth”  energy  dependence.  Several  kinetics  models  were 
investigated  of  which  the  simplest  one  to  satisfy  the  three 
conditions  is  as  follows.  At  each  energy,  phase  space  is  divided 
into  two  regions  whose  volumes  are  proportional  to  a 
parameter  F1  and  F2  =  1  —  F,.  Trajectories  in  region  1  over 
time  hop  into  region  2  with  a  rate  fc12.  Trajectories  hop  from 
region  2  to  region  1  with  a  rate  lc21  that  in  general  is  not  equal 
to  k12.  Isomerization  and  dissociation  only  occur  out  of  region  2 
with  rates  fcDIS  and  k1S0.  If  C,(f)  is  the  population  in  region  2  of 
the  ith  isomer  and  Q,(t)  the  corresponding  population  in  region 
1,  then  from  the  description  above  the  first  order  kinetics 
equations  are 

dC,-(f)/df  =  —  (^dis  "F  ^iso  +  ^21  )Q(0  +  fcisoQ— 1(0 

+  fci2Q.,(0 

dQ,.(f)/dt  =  - k12Q_.(t )  +  k2lCj(t) 

(8) 

where  CI<0(t)  =  Q,<0(f)  =  0,  C0(0)  =  Fv  and  Qo(0)  =  F,. 

Since  P,(f)  =  C,(f)  +  Q,(t),  it  is  easy  and  more  convenient  to 
convert  the  above  into  equations  in  those  for  C,(f)  and  P,(f): 

dCj(t)/dt  —  —  (^dis  "F  kiso  +  k2 1  +  kl2)Ci(t) 

+  kisoCj-i(0  +  k12fj(0 

dPj(t)/dt  =  —  (kDis  +  kiso)Q(0  +  kISoQ-i(0  (9) 

where  Ci<0(t)  =  P,<0(0  =  0,  C„(0)  =  F2,  and  P0(0)  =  1. 

The  above  equations  can  be  analytically  solved  as 


Pl(t)  =  o~at  £  A,./  +  e~ht  £  B,/ 
i= o  ;=o 

j=o  j= o  (10) 

where 

S  +  ,/i  —  4[fcDIS  +  kiso]ki2 

a  =  - 

2 

.  S  —  <J~S  —  4[fcDIS  + 

b  =  - 

2  (11) 

where  S  =  fcDIS  +  fcIS0  +  fc12  +  fc^. 

The  four  sets  of  expansion  coefficients  are  best  expressed  by 
recursion  expressions.  The  following  expressions  are  the  initial 
values  for  the  recursion: 


^oo  — 


^oo  — 


k  ^[kpis  +  kls0] 

a  —  b 

a  ~  fi2[kpis  +  kisp] 
a  —  b 
fl^oo 


^dis  +  ^-iso 


Dn 


bBn 


^-DIS  +  ^-iso 


For  i  >  0,  the  recursions  for  At-  take  the  form 


(12) 


I  n  - 


FtsoT  F  Fdis  ^isolQ-ij-i  +  kiso/Q-i,j+i  j(j  +  l)Ay+i 

j{a  -  V) 

Ai 


Al,0 


■ 

Ai  ~  ijsoF 
a  —  b 


(13) 


In  eq  13,  Aj.  |>it  =  C^i>k  =  AH<0  =  0.  The  recursion  for  B is  the 
same  form  as  eq  13  with  the  substitution  of  ( b ,  a)  for  (a,  b ),  Dkt 
for  Cy,  and  Bki  for  Akl.  For  i  >  0,  the  recursion  for  C„  is 


cr 


aAij  ~  O’  +  l)Ai,j+ 1  +  Q-l 


^-dis  "F  ^-iso 


(14) 


where  Akl>k  =  Cki>k  =  0.  The  recursion  for  D ^  is  the  same  form 
as  eq  14  with  the  substitution  of  a  for  b,  Dkl  for  Ckl,  and  Bkl  for 

Akl- 

This  kinetics  model  has  the  isomerization  and  dissociation 
processes  explicitly  represented  through  specific  rate  constants. 
The  form  of  the  analytic  solution  has  all  isomer  populations 
controlled  by  two  exponential  terms  whose  relative  weights  and 
differences  in  arguments  are  controlled  by  adjustable 
parameters.  The  model  heuristically  represents  the  common 
explanation  for  multiexponential  unimolecular  dissociation, 
namely,  that  IVR  is  too  slow  to  keep  phase  space  uniformly 
occupied  during  fast  bond-breaking  processes,  leading  to  a 
breakdown  in  the  statistical  model  responsible  for  single 
exponential  decay.  In  the  above  kinetics  model,  the  complex 
IVR  process  has  been  drastically  simplified  to  the  coupling  of 
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two  regions  of  phase  space.  Strong  coupling  allows  all  of  the 
phase  space  to  readily  participate  in  dissociation  or  isomer¬ 
ization.  Weak  coupling  means  that  the  F,  volume  of  phase 
space  lags  in  its  participation  that  in  turn  affects  the  rate  of 
dissociation  and  isomerization. 

The  kinetics  model  just  described  has  five  adjustable 
parameters,  kDIS,  klso,  kn,  k21,  and  Fj.  This  model  does 
represent  the  trajectory  data  quite  well,  but  a  simpler  model  fits 
the  data  almost  as  well.  The  simplicity  is  achieved  by  relating 
k12  and  k21  to  F1.  k12  and  k2l  can  be  thought  of  as  backward  and 
forward  rates  that  control  the  approach  to  equilibrium  between 
the  two  divided  regions  of  phase  space.  Their  ratio  can  be 
thought  of  as  an  equilibrium  constant,  and  such  constants  are 
always  ratios  of  the  phase  space  volume  of  the  two  species  that 
are  in  equilibrium.  To  the  degree  that  the  phase  space  volume 
of  regions  1  and  2  is  proportional  to  F1  and  1  —  Fv  then 


^21 


!' 


This  relationship  reduces  the  number  of  adjustable  parameters 
to  four  and  is  the  version  of  the  kinetics  model  used  to  fit  the 
trajectory  results  at  all  six  energies. 

The  four  parameters  in  the  kinetics  model  were  fit  by 
nonlinear  least-squares  techniques  to  the  trajectory  values  for 
Pi(t)  for  i  =  0,  1,  2,  3.  For  each  value  of  energy,  there  are  a  little 
more  than  2000  values  of  t  for  which  F,(t)  was  recorded  for 
each  value  of  i.  The  resulting  ~8000  data  points  were  treated 
equally  in  the  least-squares  process  to  arrive  at  optimal 
parameters  that  minimized  the  rms  error.  Initial  runs  showed 
that  the  optimal  value  of  F1  varied  significantly  from  energy  to 
energy,  producing  optimal  values  in  the  other  parameters  that 
had  an  irregular  dependence  on  energy.  However,  near  optimal 
values  of  F1  insignificantly  raised  the  rms  error.  Imposing  on  F1 
various  low  order  polynomial  dependencies  in  E  indicated  that 
a  constant  F1  with  respect  to  E  was  both  the  simplest  and 
optimal  representation.  The  final  least-squares  solution 
optimized  kDIS,  fcIS0,  and  fc12  independently  at  each  energy  to 
minimize  the  rms  error  at  that  energy  but  optimized  a  constant 
value  of  Fj  to  minimize  the  six  term  sum  of  the  rms  errors  at 
each  energy. 

In  Figure  4,  the  dash— dot  lines  represent  the  resulting 
optimal  fit  for  three  of  the  six  energies.  The  other  three  energies 
are  quite  comparable.  The  rms  errors  range  from  0.0014  (E  = 
120  kcal/mol)  to  0.0028  (F  =  150  kcal/mol).  Since  the 
populations  are  all  normalized,  this  scale  of  rms  errors  indicated 
overall  fits  at  the  few  tenths  of  a  percent  level.  Figure  4 
indicates  that  the  emergence  of  multiexponential  decay  at  the 
higher  energies  in  P0(t )  is  captured  well  by  this  model. 
Although  the  kinetics  model  results  are  not  plotted  in  Figure  5 
that  highlights  the  emergence  of  multiexponential  behavior, 
they  would  be  nearly  indistinguishable  from  the  trajectory 
results.  The  energy  dependence  of  the  optimal  values  of  fcDIS, 
fcIS0,  and  fc12  is  displayed  by  the  symbols  in  Figure  9;  the 
optimal  value  of  F1  is  0.71.  The  energy  behavior  of  all  three 
parameters  is  reasonably  smooth  and  regular.  The  solid  lines  in 
Figure  9  are  a  least-squares  fit  of  the  natural  log  of  the 
parameters  to  a  quadratic  in  the  energy.  fcDIS  and  k1S0  are 
exceptionally  well  represented,  but  there  are  some  relative 
minor  deviations  in  the  fit  of  fc12. 

The  results  in  Figure  9  confirm  the  expected  behavior.  At  low 
energies,  k12  kDIS  +  kIS0  and  a  single  exponential  dominates 
the  decay.  Analysis  of  eq  11  shows  that  the  exponential 


E  (kcal/mol) 


Figure  9.  Least-squares  values  of  the  three  rate  constants  in  the 
kinetics  model  versus  energy  (symbols).  Least-squares  fits  of  the 
energy  dependence  of  the  optimal  values  of  the  three  rate  constants 
(solid  lines). 


arguments  a  and  b  are  proportional  to  k12  and  fcDIS  +  fcIS0, 
respectively,  but  according  to  eq  13,  the  faster  decay  receives 
very  little  amplitude  and  is  just  a  short-time  transient.  At  high 
energies,  ku,  fcDIS,  and  fcIS0  are  all  comparable  and  a 
biexponential  decay  emerges.  From  the  trends  in  Figure  9,  at 
higher  energies,  still  kl2  <SC  fcDIS  +  fcIS0.  Then,  a  and  b  will  retain 
the  same  proportionality  to  the  rates  but  the  amplitudes  of 
those  terms  will  be  proportional  to  F2  and  Fv  respectively.  For 
the  optimal  Fj  value  of  0.71,  biexponential  unimolecular  decay 
would  be  quite  evident.  However,  nothing  in  this  analysis  of  our 
trajectory  results  points  to  an  extremely  slow  rate.  Rather,  all 
rates  increase  with  E  but  not  equally  fast.  The  slow  rates 
measured  by  the  Chen  and  Fischer  groups  would  require  rates 
orders  of  magnitude  lower  than  those  in  Figure  9. 

The  results  in  Figure  9  indicate  that  fcDIS  and  fcIS0  are 
comparable  throughout  the  energy  range  of  the  trajectories. 
Considering  that  the  barriers  to  dissociation  and  isomerization 
are  quite  close  to  one  another  (relative  to  values  of  E  more  than 
twice  either  barrier),  this  is  perhaps  not  a  surprising  result. 
When  examined  more  quantitatively,  the  small  differences 
between  fcDIS  and  fcIS0  do  have  consequences.  As  seen  in  Figure 
9,  there  is  a  finite  energy  range  over  which  fcIS0  >  fcnis-  As 
discussed  earlier,  the  more  isomerization  dominates  dissocia¬ 
tion,  the  lower  FH  becomes.  Because  the  relative  size  of  fcIS0  to 
kDIS  is  not  monotonic  in  energy,  Fh  will  also  not  be  monotonic 
in  energy.  This  is  illustrated  by  the  black  line  in  Figure  7  where 
Fh  is  calculated  over  an  extended  energy  range  using  the  energy 
dependence  of  the  fit  in  Figure  9  and  the  statistical  alternation 
in  Figure  8.  The  minimum  in  Fh  corresponds  to  the  maximum 
in  fcIS0  relative  to  kDIS  in  Figure  9.  The  blue  triangles  in  Figure  7 
involve  only  an  interpolation  of  trajectory-calculated  alterna¬ 
tions  in  Figure  8.  They  also  display  a  shallow  minimum, 
indicating  that  the  fit  is  representing  a  real  feature  of  the 
trajectory  results.  This  suggests  that  precise  measurements  of 
the  energy  dependence  of  isotopic  scrambling  could  reveal 
subtle  variations  in  the  competition  between  isomerization  and 
dissociation. 

The  energy  dependence  of  the  rates  as  displayed  in  Figure  9 
has  a  limited  range  of  validity  when  extrapolated  beyond  the 
energy  range  of  the  trajectory  calculations.  At  the  low  energy 
end,  the  higher  barrier  to  isomerization  relative  to  dissociation 
means  that  isomerization  rates  should  become  vanishingly  small 
as  E  approaches  the  isomerization  barrier  height  of  ~43  kcal/ 
mol  (see  Table  l).  The  fitted  result  is  small  but  not  small 
enough  to  result  in  Fh  equaling  the  required  value  of  1,  as  seen 
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in  Figure  7.  At  the  high-energy  end,  a  few  tens  of  kcal/mol 
beyond  the  energy  range  of  Figure  9,  fcIS0  begins  to  decrease 
with  energy,  a  physically  unlikely  scenario. 

Within  the  energy  range  of  the  trajectories,  the  fit  is  a  close 
representation  of  the  trajectory  results  as  seen  in  Figure  4  but 
underestimates  Fh  as  seen  by  the  black  curve  in  Figure  9  by 
<0.05.  The  majority  of  that  underestimation  is  due  to  the  use 
of  statistical  alternations  that  are  shown  in  Figure  8  to 
underestimate  trajectory  results.  However,  at  the  higher 
energies,  there  is  a  residual  underestimation  of  up  to  0.02 
that  comes  from  an  imperfect  fit.  As  indicated  at  the  beginning 
of  this  section,  more  elaborate  kinetics  models  are  possible  that 
would  further  reduce  fitting  error.  However,  it  is  important  to 
note  that  the  fitting  target  was  the  evolution  of  the  C2H5 
population,  not  the  H  product  population,  with  time.  One 
could  have  fit  the  time  dependent  population  of  H  from 
dissociation  after  each  isomerization.  Since  is  all  about 
products,  it  is  likely  such  a  fitting  strategy  would  have  improved 
agreement  with  trajectory  calculations  of  Fh  at  some  modest 
cost  to  representing  trajectory  populations  after  each  isomer¬ 
ization. 

5.  CONCLUSIONS 

Motivated  by  photodissociation  experiments  on  the  ethyl 
radical  that  have  discovered  very  slow  dissociation  processes, 
we  have  performed  a  trajectory  study  of  the  dissociation  and 
isomerization  of  C2HS  out  to  as  much  as  100  ps  over  an  energy 
range  of  100—150  kcal/mol.  The  feasible  propagation  of 
trajectories  for  such  a  long  time  requires  an  efficient  model  for 
the  potential  energy  surface.  In  the  first  gas-phase  application 
we  are  aware  of,  we  used  the  AIREBO  empirical  molecular 
mechanics  potential.  In  its  native  form,  AIREBO  has  both 
technical  and  accuracy  limitation  in  representing.  However,  a 
limited  generalization  of  the  potential  energy  form  and  a 
limited  parameter  optimization  to  select  calculated  CBS-QB3 
saddle  point  properties  for  both  dissociation  and  isomerization 
leads  to  a  realistic  AIREBO  specific  reaction  parameter  (SRP) 
surface  appropriate  for  the  unimolecular  decay  of  C2H5.  The 
resulting  AIREBO-SRP  surface  largely  retains  the  computa¬ 
tional  speed  of  the  native  AIREBO  surfaces,  allowing  the 
feasible  calculation  of  very  long  time  trajectories. 

At  each  energy,  10  000  trajectories  were  followed  until 
essentially  all  dissociated.  This  is  an  order  of  magnitude  longer 
in  time  and  larger  in  sample  size  than  any  previous  study.  The 
trajectory  results  show  no  evidence  of  a  very  slow  dissociation 
process.  They  do  show  evidence  of  multiexponential  decay  at 
the  higher  end  of  the  energy  range,  but  the  rate  constant 
components  of  this  decay  are  still  many  orders  of  magnitude 
faster  than  slow  decay  observed  and  attributed  to  unimolecular 
processes  on  the  ground  state  surface.  The  fraction  of  product 
energy  deposited  in  translation  is  about  75%  of  the  most  recent 
measured  value.  Isotopic  scrambling  is  overestimated  for 
CD2CH3,  with  the  fraction  of  products  constituting  H  + 
CD2H2  calculated  at  about  ~0.8,  while  the  measured  value  is 
0.95  ±  0.05.  The  alternation  with  repeated  isomerization  of  the 
H  fraction  of  products  can  be  largely,  but  not  completely, 
explained  by  an  isomerization  model. 

A  divided  phase  space  kinetics  model  of  dissociation  and 
isomerization  fits  the  trajectory  results  quite  well.  The  divided 
phase  space  aspect  of  the  model  has  one  region  of  phase  space 
that  can  isomerize/dissociate  and  a  second  region  cannot.  The 
IVR-like  rate  of  transfer  between  those  two  regions  relative  to 
the  dissociation  and  isomerization  rate  can  influence  product 


formation.  The  model  represents  the  growth  of  multi¬ 
exponential  decay  with  increased  energy,  as  seen  in  the 
trajectory  results.  The  model  represents  in  detail  repeated 
isomerization  processes.  Interestingly,  the  model  has  isomer¬ 
ization  and  dissociation  rates  that  vary  relative  to  one  another 
in  a  non-monotonic  fashion  as  a  function  of  energy,  resulting  in 
a  H  fraction  of  products  that  has  a  minimum  with  respect  to 
energy.  This  suggests  that  precise  measurements  of  the  energy 
dependence  of  isotopic  product  fractions  may  reveal  subtle 
aspects  of  the  competition  between  isomerization  and 
dissociation. 

The  AIREBO-SRP  potential  in  this  work  bring  to  four  the 
number  of  potential  energy  surfaces  that  have  been  used  in 
trajectory  studies  relevant  to  photodissociation  experiments. 
Only  one  of  those  potential  energy  surfaces25  has  produced 
“trapped  trajectories”  that  indicate  a  slower  process  for  ethyl 
radical  dissociation.  However,  the  actual  value  of  the  slower  rate 
could  not  be  determined  because  of  the  large  sample  error  of 
the  handful  of  trajectories  involved  and  because  of  the  very 
short  times  (~2  ps)  the  trajectories  were  followed.  These 
results  could  be  consistent  with  our  trajectory  results  at  higher 
energies  that  show  multiexponential  decay  but  with  rates 
involved  that  are  orders  of  magnitude  different  from  the 
experimental  observations.  The  authors  of  that  study  have 
noted  the  sensitivity  of  the  trapping  to  details  of  the  potential 
energy  surface.27  Every  one  of  the  four  surfaces  has  been  a 
compromise  to  higher-level  descriptions  in  order  to  be 
computationally  inexpensive  enough  to  enable  trajectory  runs 
to  the  time  scales  necessary  to  make  contact  with  experiment. 
From  a  theoretical  perspective,  any  further  dynamics  studies 
would  best  be  performed  on  a  high  fidelity  fit  to  a  high-level 
electronic  structure  description  of  the  ground  state  potential 
energy  surface.  With  15  internal  degrees  of  freedom,  this  is  a 
challenging  task,  especially  up  to  the  energies  sampled  by  the 
photodissociation  experiments.  We  have  been  interested  in  the 
interpolated  moving  least-square  (IMLS)  method45  of  fitting 
potential  energy  surfaces  automatically  with  a  minimization  of 
the  number  of  electronics  structure  calculations  that  can  be 
performed.  Fitting  the  C2H5  surface  with  this  method  or  other 
methods  capable  of  automatic  potential  energy  generation  is  a 
direction  we  intend  to  pursue  for  future  research. 
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■  GLOSSARY 

AIREBO:  adaptive  intermole cular  reactive  empirical  bond- 
order 

CBS:  complete  basis  set 
CBS-QB3:  complete  basis  set-quadratic  Becke3 
CCSD(T):  coupled  cluster  with  single,  double,  and 
perturbative  triple 

C,(f):  population  in  region  2  of  the  ith  isomer 

DFT :  density  functional  theory 

iq:  phase  space  volume  for  region  1 

F2:  phase  space  volume  for  region  2,  F2  =  1  —  F1 

Fh-  fraction  H  atom  dissociation 

{ fT ):  fraction  of  product  energy  released  into  translation 
HCTH147(3)6-31+G**/6-31G**:  HCTH147(a)6-31+G** 
density  functional  parametrized  by  Boese  et  al.44  with  the 
6-3 1G**  basis  set 

IMLS:  interpolating  moving  least  square 

IVR:  intramolecular  vibrational  redistribution 

kl2:  rate  of  trajectories  hop  from  region  1  to  region  2 

k2l:  rate  of  trajectories  hop  from  region  2  to  region  1 

kDls:  dissociation  rate 

kiSO :  isomerization  rate 

LJ:  Lennard-Jones 

M05-2X/CBSB7:  M05-2X  density  functional  described  in  ref 
33  with  the  CBSB7  basis  set;  the  CBSB7  has  the  form  6- 
311G(2d,d,p)  used  by  the  CBS-QB3  method45 
MR-CISD:  multireference  configuration  interaction  with 
singles  and  doubles 
PES:  potential  energy  surface 
P,(f):  time-dependent  normalized  population 
PM6:  parametric  method  number  6 
POL- Cl:  polarization  configuration  interaction 
QCISD(T):  quadratic  configuration  interaction  with  single, 
double,  and  perturbative  triple 
Qj(t):  population  in  region  1  of  the  ith  isomer 
REBO:  reactive  empirical  bond-order 
RMS:  root  mean  square 
RRKM:  Rice— Ramsperger— Kassel— Marcus 
SA-MCSCF:  state-averaged  multiconfiguration  self-consis¬ 
tent  field 

SEMO:  semi-empirical  molecular  orbital 
SRP:  specific  reaction  parameters 
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